function [ U, V, err ] = pmf( Y, Omega, K, mu, lambda)
% Omega: observed entries
% |Omega(Y - U'*V)|_F^2 + mu*|U|_F^2 + lambda*|V|_F^2

M = size(Y, 1);
N = size(Y, 2);

U = randn(K, M);
V = randn(K, N);
maxIter = 10000;
stepSize = svds(Y + U'*V, 1); stepSize = 1/stepSize;
err = zeros(maxIter, 1);
% gradient descent
for i = 1:maxIter
    gradU = V*((U'*V - Y).*Omega)' + mu*U;
    U = U - stepSize*gradU;
    
    gradV = U*(Omega.*(U'*V - Y)) + lambda*V;
    V = V - stepSize*gradV;
    
    err(i) = norm((Y - U'*V).*Omega, 'inf');
    if(i > 10 && abs(err(i) - err(i-1)) < 1e-6)
        break;
    end
end
err = err(1:i);

end

